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Abstract 

We present a simple model of a dynamical system driven by externally- 
imposed coherent noise. Although the system never becomes critical in the 
sense of possessing spatial correlations of arbitrarily long range, it does orga- 
nize into a stationary state characterized by avalanches with a power-law size 
distribution. We explain the behavior of the model within a time-averaged 
approximation, and discuss its potential connection to the dynamics of earth- 
quakes, the Gutenberg-Richter law, and to recent experiments on avalanches 
in rice piles. 
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I. INTRODUCTION 



There has in the last few years been considerable interest in extended systems which 
self-organize into a state exhibiting large scale fluctuations and intermittent dynamics. One 
of the earliest attempts to model systems of this type was made in 1987 by Bak, Tang, and 
Wiesenfeld, who proposed a simple lattice model for the avalanches produced by depositing 
grains of sand on an ever-growing sand pile PJ. Despite having only short-range interactions 
and no tunable parameters, their model organizes itself into a state with long-range spatial 
correlations and avalanches of size not limited by any finite correlation length. It has been 
proposed that similar self-organized critical (SOC) behavior could lie behind a wide range of 
physical phenomena showing 1/f noise and scale-free fluctuation distributions. SOC models 
have been put forward to describe the dynamics of earthquakes , biological evolution || 
and extinction |4]], forest fires ||, and many other systems ||. The common features of these 
models are that (i) they are all driven very slowly, and (ii) they all have perfect memory, 
i.e., in the absence of the driving force the model would be entirely stationary. 

A distinctive observable consequence of SOC dynamics is that the distribution of fluctu- 
ation (or avalanche) sizes takes a power-law form with characteristic exponent r: 

PavalO) OC S~ T . (1) 

The value of r typically lies in the range 1 < r < |, with the value | corresponding to a 
critical branching process, appearing if one makes the "random neighbor approximation" 
in which each site interacts with a randomly-selected small number of other sites 0. This 
approximation is equivalent to the limit of infinite dimension, and should give correct results 
for systems above their critical dimension. In reality however, a number of the systems which 
are modeled using SOC dynamics in fact display event size distributions with fairly large 
exponents. Terrestrial earthquakes, for example, appear to follow the Gutenberg-Richter 



law || with r « 2.0, and the ID rice pile experiment of Frette and co-workers ||10|| , which 



has been compared to the sandpile model of Bak, Tang and Wiesenfeld, gives r = 2.1 ± 0.1. 
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A number of models have been proposed which offer explanations for these higher values of 
t. The model discussed by Christensen and Olami is one such, but it achieves its result 
only at the price of having an entirely deterministic dynamics; if one introduces randomness 
into the model, the simple scaling behavior is destroyed. The reason is that when the 
exponent describing the distribution of avalanches' spatial extent becomes larger than 2, the 
mean avalanche size becomes finite and independent of system size, and the spatial overlap 
between subsequent avalanches becomes insignificant |§ . In the presence of randomness, this 
can prevent the system from building up any long-range correlations, and ultimately destroy 
the critical state. We conjecture that, in this regime, any randomness in the positions of 
the nucleation centers of the avalanches will destroy self-organization of long range spatial 
correlations. 

In this paper, we present a different explanation to account for systems that have larger 
values of r. We demonstrate that power-law event size distributions having r around 2 or 
greater, are typical of extended systems with quenched memory if they are driven by coherent 
noise, and that in such systems they are present even in the absence of any interaction 
between the different parts of the system. (This is different from the situation in the SOC 
models, where the system is driven by a local driving force, coupled with interactions between 
the components of the system.) The simplest model demonstrating the phenomenon is 
defined as follows. Consider a system of N agents, such as grains on the surface of a sand 
pile or points of contact in a subterranean fault. With each agent i we associate a threshold 
for movement Xj which can take values falling in some specified range and represents the 
amount of stress that the agent will withstand before it moves. For convenience, we choose 
to measure Xi on a scale on which < Xi < 1. The dynamics of the model then consists of 
the repetition of two steps: 

1. A fixed fraction / of the agents are selected at random, and the values of their threshold 
variables are exchanged with new random numbers selected uniformly from the 
interval [0, 1). 
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2. We select a random number or "stress level" 77 from some distribution Pstress(^)- All 
Xi below f) are also exchanged with new random numbers selected uniformly from the 
interval [0, 1). The number of agents whose thresholds are changed in this fashion is 
the size s of the avalanche taking place in this time-step. 

The random selection of different values for 77 at each step may be thought of as imposing 
external stresses which coherently (in other words, simultaneously) influence all of the weaker 
agents — those having suitably low thresholds for stress — but leaves unchanged the stronger 
ones. It seems physically reasonable to assume that smaller stresses should be more common 
than larger ones, and in the following discussion we make the assumption that p s tress(?7) is 
largest at 77 = and falls off to zero as 77 becomes large. We denote the typical scale of the 
falloff by a. The most interesting regime is when o <C 1 and / -C 1, 

II. RESULTS 

We have examined the properties of this model both analytically and numerically. In- 
stead of simulating the model directly, we have developed an algorithm which calculates 
the threshold distribution and avalanche sizes in a formally exact way for a system with 
N = 00. Starting off with a uniform distribution of thresholds, the system evolves towards 
a statistically stationary state. In this state we record the mean threshold distribution and 
the frequency distribution of avalanches. The results are shown in Figures [l] and ^| for a 
simulation using exponentially distributed stresses: 

Pstressiv) = -eXp(-77/cr). (2) 

a 

As Figure [I] shows, the distribution of avalanche sizes s is flat up to a certain point (whose 
position varies with a and /) and then falls off as a power law according to Eq. ([[]) with 
t Pd 2.0. This power-law behavior appears to be robust in the regime of small / and a. If, 
for example, instead of Eq. (^j we employ a Gaussian stress distribution then, although the 
average distribution of thresholds (Figure |2|) changes radically, the power-law form of the 
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avalanche distribution remains. Notice, however, that the exponent r changes slightly as the 
applied stresses are varied. For the Gaussian distribution, for example, we find r = 2.2 ±0.1, 
as opposed to r = 1.9 ± 0.1 for the exponential. And for steeper distributions of stresses 
(p(rf) oc exp{—{rj/a) q ) with q > 4) we find r = 2.4 or greater. 

In order to investigate possible connections with spatially-organized models, we have also 
implemented our model on a lattice and at each time-step eliminated not only those agents 
whose thresholds for stress fall below the selected level, but also their neighbors. In all cases 
we observe a power- law distribution of avalanches with exponent in the vicinity of r = 2. 

In order to understand the appearance of this power law, let us consider the time- 
averaged behavior of the model. The statistically stationary state arises as a competition 
between the two processes comprising the dynamics: the stresses which tend to remove 
lower thresholds from the distribution and thus shift the weight of the distribution to higher 
values of x, and the aging, which tends to move weight back down again. The result is 
that the average threshold distribution pthresh(^) is a highly nonhomogeneous, monotonic 
increasing function of x which, for small a, tends to have a plateau as x approaches unity 
(see Figure |2|). By balancing the two competing processes, we can calculate Pthresh^) and 
hence the avalanche distribution. For concreteness, we perform the calculation here for the 
exponentially distributed stresses of Eq. @. 

The probability of an agent possessing a threshold x lying below the stress level r\ at any 
given time-step (and hence of it moving during this time-step) is 



The total time-averaged rate at which agents move in the interval between x and x + dx is 
then 



where the ^-independent constant W on the right-hand side is the time-averaged rate at 
which probability is added to pthresh- Rearranging we have 




(3) 



Pmove{%) Pthresh(^ ) dx + /pthresh(z) dx = W dx 



(4) 
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The constant is easily fixed by requiring that pthresh(^) integrate to unity giving 



(5) 
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a 
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Je 1 ^ + 1 

7+i 



(6) 



For small / and a, pthresh(^) rises exponentially from zero and then levels off in a plateau 
around x = —a log /. Physically this arises because agents possessing thresholds above this 
point are affected only by the aging process, which treats them all equally. Below this level, 
the stress process is important too, and it preferentially moves those with lower thresholds. 
The avalanche size distribution is given by 



P(S\V) PstressO?) cfy. 



(?) 



The probability p{s\r]) of getting an avalanche of a certain size given a certain stress level, 
depends on the distribution of thresholds, which will in general vary from one time-step 
to another. However, if we make the "time averaged approximation" (TAA) whereby one 
assumes that at each time-step the threshold distribution can be approximated by its time- 
averaged value, then p(s\r)) = 6(s(r)) — s) where 5(77) is just 

™ 

s (v) = / PthrcshO) dx. (8) 

JO 

The avalanche size distribution then becomes 



PavallSj 



8(s(rj) - S) Pstress(^) dl] = PsbreBB ^^ 
PthreshlMs)) 
1 
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(/ + e-"W/ ff ) 



(9) 



where we have used Eqs. © and (§). We can calculate the stress level r)(s) corresponding 
to an avalanche of size s from the same two equations, which give 



1 + f e ^ a 
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l + f 



(10) 



6 



for e _1 /°" <C / 1 and a C 1. We can now distinguish a number of different regimes. For 
small avalanches, such that s <C a, the logarithm on the right-hand side can be expanded 
giving s + of 7a afe v ^ a . Substituting into Eq. (§) 

Paval(s) OC [S + af]~ 2 for S < (7. (f f ) 

This gives a flat avalanche distribution for small s up to about s = af, and then a power- 
law distribution for larger s with exponent r = 2. The approximation breaks down when 
s ~ a, giving way to a regime in which e n ^ a ~ e s , and hence the avalanche distribution falls 
off exponentially with s. The various regimes can clearly be seen in the numerical results 
presented in Figure [l], and the predicted cross-over points between them agree well with the 
theory. 

When / decreases below e" 1 / ", the approximations in Eq. fllOl) break down and instead 
it becomes valid to write e v ^ a wl + se 1 ^. In this regime the theory predicts a breakdown in 
the scaling, a phenomenon which is also seen in the simulations. Thus the reloading process, 
whose scale is set by /, must be small but necessarily non-zero if we are to see power-law 
behavior in the avalanche distribution. Notice however that at precisely / = the theory 
predicts a return to r = 2 scaling, which is not seen in the simulations, implying that the 
TAA breaks down in this regime because the distribution p(s\rj) becomes too broad to be 
well approximated by a 5-function. 

The physical principle behind the appearance of a power-law distribution here is the 
interdependence of the avalanche and threshold distributions; the avalanche distribution 
is a function of the particular distribution of thresholds at any time, but the threshold 
distribution is itself produced by the action of the avalanches. 

III. CONNECTION WITH OTHER MODELS 

There are clear similarities between our model and the sand pile model, in which sites also 
possess a certain threshold stress that they will withstand without adjusting. Furthermore 
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we have a source term, the reloading or aging fraction / of agents which at each time-step 
loose memory of their previously assigned thresholds. This source term is similar in effect 
to the addition of the single grains of sand in the sand pile models. There are however 
some important differences between our model and the SOC models. First, the stresses 
in our model are coherent, rather than localized as they are in the sandpile. Second, the 
agents are, at least in the simplest versions of the model, entirely non-interacting. In SOC 
models, it is the interactions which give rise to avalanches. In our model on the other 
hand the avalanches of simultaneously moving agents arise because all the agents feel the 
same externally imposed stresses. There is no causal connection between the events which 
comprise an avalanche; each agent moves independently of the others. 

Unlike other model systems for large scale fluctuations, such as the Burridge- 



Knopoff (BK) model JTTJ and the recycled version of the Democratic Fiber Bundle 
Model (DFBM) fL2| , the model presented here does not make a clear distinction between 
small, finite-sized events, and large ones whose size scales like the size of the system. In 
the BK model, for instance, the spectrum of event sizes contains two separate parts, one 
composed of small events which scales as s~ 2 , and another composed of the big events, which 
occur quasiperiodically. The BK and DFBM models are not statistically stationary, by con- 
trast with our model whose dynamics rapidly reaches a statistically stationary state. Models 
such as BK and DFBM also show "foreshock" events in which large avalanches are preceded 
by smaller ones. Our dynamics does not have foreshocks but does display aftershock events, 
a phenomenon which we discuss in greater detail in next section. 



IV. DISCUSSION 

Next, we would like to examine the potential relationship of our model to processes 
occurring in real physical systems. First we consider earthquakes. To begin with, we ignore 
spatial correlations and consider the variables xi to be thresholds for movement at various 
points along a fault. The coherent stress r\ is provided by long- wavelength background noise 
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from some external source, such as other distant tremors, or movements in the deeper regions 
of the earth, and the reloading / is due to slow plastic deformation from tectonic movements 
of the crust. As we have seen, these elements alone lead directly to a power-law distribution 
of earthquake sizes very close to the observed Gutenberg-Richter law, without the need to 
invoke interactions between neighboring parts of the fault. That is not to say that such 
interactions do not exist, only that they are not necessary to produce the observed power 
law. (Kagan |13| has presented evidence of a fractal pattern in the spatial distribution of 



earthquake activity, which is an indication that interactions are a feature of the dynamics. 
This however need not lead us to conclude that these local interactions are necessary for 
producing the observed size distribution of events.) 

Another interesting feature of our model is that it shows clear aftershocks. The mecha- 
nism for these is straightforward. When a large avalanche takes place, a significant fraction 
of the thresholds in the system are replaced with new, uniformly distributed ones. Be- 
cause of the monotonic increasing form of the threshold distribution, this has the effect of 
shifting the weight of the distribution downwards, increasing the fraction of agents with 
low threshold for movement. The result is that subsequent stresses on the system have a 
larger-than-normal effect, and we see an amplification of the usual level of "background" 
avalanches in the aftermath of a particularly large event. In Figure we show a section 
of a time series of avalanches from one of our simulations, which clearly displays this af- 
tershock effect. Notice that if we apply the argument iteratively, we would also expect to 
see sequences of "after- aftershocks" following each of the aftershocks, a behavior which is 
indeed evident in Figure [| We have also measured the average probability of getting an 
event of significant size in the aftermath of another large one, and found that for small times 
its distribution goes approximately as t _1 (Figure ^). A similar result is seen in the data 



from real earthquakes, and is commonly referred to as Omori's law |L4|. 

The t~ x distribution can be understood as follows. A large avalanche will redistribute 
the thresholds of a large fraction of the agents in the system uniformly across the interval 
of allowed values (0 < x < 1 in this case). A subsequent stress of magnitude rji will remove 



all those agents with x < rji, and produce an aftershock extinction of a certain magnitude. 
In order to get another significant aftershock we now need a stress 7] 2 > r]i in order to reach 
those agents which were not affected by the first aftershock. In general, if it took a time t 
to get the first stress, then it will on average, take as long again to get another of the same 
magnitude, or an aggregate time of 2t until a stress of size 772 comes along. Repeating the 
argument, it will take as long again, or a total time of At to get the third aftershock, and 
so forth. Given this exponential increase in the time intervals between these events, it is 
not hard to show that the histogram of aftershock events should have a t _1 power-law form, 
regardless of the precise distribution of stresses applied to the system. 

Note that our mechanism is by no means the only way to obtain aftershocks. An alterna- 
tive mechanism has been proposed by Nakanishi [[OJ using a Burridge-Knopoff-like model in 



which relaxation processes are introduced by considering the geometry of stress redistribu- 
tion following large quakes. As with the BK model, Nakanishi's model has a quasiperiodic 
dynamics. 

Second, let us compare our model with the results of recent studies of one-dimensional rice 



piles by Frette et al. |TD| . In these studies the experimenters found a frequency distribution 
of avalanche sizes s which was flat up to a certain fraction of the total size of the pile, and 
then fell off as a power of s for larger avalanches according to Eq. ([!]), with a measured 
exponent of r = 2.1 ± 0.1. A similar behavior is seen in the simulation results from our 
model (Figure |l|) which also display a flat distribution of avalanches up to a certain fraction 
~ erf of the total system size, and then a power-law fall in avalanche frequency for larger 
sizes with exponent close to two. A possible interpretation of the experimental data then is 
that the dynamics of the rice pile is one of avalanches produced by the interplay of reloading 
with coherent stress. The reloading / could arise as a result of newly added grains of rice, 
which tend to randomize the thresholds for grains on the surface, and the stresses might 
come from the tumbling of new grains as they are added to the pile. The plateau in the 
avalanche distribution for small sizes s is then caused by rice grains which tumble past a 
number of sites before coming to rest, but have only enough energy to disturb the most 
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unstable of those sites, and the larger events which form the r = 2 power law are the result 
of occasional larger stresses in the tail of the distribution. The one-dimensional nature of the 
system ensures that all input disturbances propagate through a large portion of the system, 
and thus may be treated as coherent. In a two-dimensional system this would not be the 
case, and the pile might well show entirely different dynamics, either possessing a shallower 
power-law distribution r < 2, indicating perhaps that a true SOC dynamics is at work, or 
not possessing a power-law distribution at all, indicating that coherent driving forces are 
the only mechanism responsible for power laws in this system. 

Our model also makes quantitative predictions about the scaling of the line between 
the two regimes in the avalanche distribution: the position of the line should go like Naf, 
the factor of the system size N appearing when we shift from measuring avalanches as 
fractions of the system size to measuring the total energy they release. Scaling of precisely 
this form with iV is indeed seen in the experiments. Frette et al. also mention that simple 
scaling disappears when the experiment is repeated with "rounder" rice. We can explain 
this result in terms of the narrower distribution of thresholds that round rice can support, 
which corresponds to larger values of both a and /. 

The results of these experiments have also been modeled by Christensen [HJ using a 



SOC model with interacting elements. Clearly, there are aspects of the dynamics captured 
by their model which are missing from ours, particularly geometrical effects concerned with 
the spatial distributions of avalanches and the corresponding transport properties of rice in 
the pile [|17j]. However, because the exponent r is greater than 2, making (s) independent 
of system size, we can expect these properties to be independent of the largest avalanche 
events (though on the other hand, they should now depend strongly on the position and 
nature of the crossover between the two regimes of the avalanche distribution). We suggest 
that the reverse is also true, i.e., that the observed large avalanches could appear even in 
the absence of long range spatial correlations. 

One characteristic which does seem to distinguish our model from the SOC alternatives 
is the existence of aftershock events. It might therefore might be profitable to investigate the 
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existence of aftershock avalanches in the experimental data, in order to make a quantitative 
distinction between the two classes of dynamics. 

Finally we would like to point out that models of the type introduced here do not 
constrain r to values close to 2. Although the values found with the simple version of the 
model outlined in Section II all lie approximately in the range 1.8 < r < 2.4, we have 
investigated other variants on the model which produce values outside this range. One 
particularly interesting version is one in which we allow for the possibility of there being 
many different kinds of stress on an agent. We suppose that agent i is subject to M 
independent types of stress, and that it has a separate threshold for yielding to each one, 
making Xj an M-dimensional vector quantity. One then assumes that all M components 
of Xj are to be replaced with new values every time any one of the types of stress exceeds 
the corresponding threshold value. In the limit M = 1, this model is clearly just the same 
as the version discussed above, and for higher values of M we continue to see a power-law 
distribution of avalanche sizes, regardless of the nature of the applied stresses. However, 
the exponent of the power law becomes steeper as the value of M increases, and appears 
to approach 3 as M becomes large. (We have investigated the model numerically up to 
M = 50.) It is interesting to note that stock market fluctuations show power-law fluctuation 
distributions with exponents close to r = 3 |1| . One may speculate whether these so-called 
"fat tails" in the distribution are the natural response to the action of external stresses on 
the market (of which there are indeed many). 



V. CONCLUSION 

To summarize, we have demonstrated that coherent noise in large systems typically gives 
rise to intermittent behavior with an "avalanche" type dynamics characterized by a power- 
law distribution of avalanche sizes with exponent in the vicinity of r = 2. This value is similar 
to that seen in a number of real systems, including rice piles and earthquakes, suggesting 
that these systems may in fact be driven by external noise, rather than self-organizing under 
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the influence of short-range internal interactions. If one allows more elaborate types of stress 
on the system one can obtain power laws with exponents as high as r = 3. 

We believe that the study of systems driven in this fashion by coherent external noise may 
offer new interpretations of intermittent dynamics in a variety of extended non-equilibrium 
systems in terms of a direct interplay between small scale structures and long wavelength 
fluctuations in the system. Such systems might include not only the ricepiles and earthquakes 
considered here, but possibly also extended chaotic systems such as economics |L8f and 



turbulence 19 
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FIG. 1. Simulation results for the frequency distribution of avalanches with exponentially dis- 
tributed stresses (solid line) and Gaussian ones (dashed line), with / = 10~ 3 and er = in each 
case. 
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FIG. 2. Simulation results for the time-averaged distribution of thresholds x with exponentially 
distributed stresses (solid line) and Gaussian ones (dashed line). As in Figure [[],/ = 10 -3 and 
a = in each case. 

16 



0.3 



tn 

CD 
N 
- 1 - 1 

(7) 



0.2 - 



<d 
o 

a 

£ 0.1 

CO 

> 

CO 



0.0 



n — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — r 



1 I I I L 



.LLa. 



_L I I L_ 



200 400 600 800 1000 

time t 

FIG. 3. Time series plot of avalanche sizes during a portion of a simulation showing clear 
aftershocks. 
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FIG. 4. Histogram of the time distribution of aftershocks following a major avalanche. The 
histogram follows a power law with an exponent close to one (Omori's law). 
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